The packages we will use are included below:

library(tidyverse) # Used for general graphing
library(ggvoronoi) # Creates Voronoi diagrams
library(deldir) # Obtains info on diagrams
library(ggmap) # Graphs onto a street map of NY
library(sp) # For heat maps

The ggvoronoi and deldir packages are what we will use to create the Voronoi diagrams. The other three packages will be used to plot the diagrams.

For our project we created numerous functions to create our diagrams and plot them.

Our first function, The function get_area will return the areas of the Voronoi polygons so that we can compare them.

get_area <- function(df){                          # Takes in data set and finds areas of 
  areas <- deldir(df$Lon,df$Lat)$summary$dir.area  # voronoi polygons for numerical analysis
  return(areas)
}

The following functions are used to easily access coordinates used to make Voronoi polygons

evparse <- function(x){         # A simple function to read a string as code
  eval(parse(text = paste0(x)))
}

get_poly_x <- function(df){                  # Obtains the x coordinates of voronoi polygons
  tiles <- deldir(df$Lon, df$Lat) %>% 
    tile.list
  len <- seq_along(tiles)
  temp_vec <- paste0("tiles$pt.", len, "$x")
  return(map(temp_vec, evparse))
}

get_poly_y <- function(df){                  # Obtains the y coordinates of voronoi polygons
  tiles <- deldir(df$Lon, df$Lat) %>% 
    tile.list
  len <- seq_along(tiles)
  temp_vec <- paste0("tiles$pt.", len, "$y")
  return(map(temp_vec, evparse))
}

The get_coords and gen_map are both designed to obtain a map or New York city centered on the data

get_coords <- function(x_vals,y_vals){     #Takes in the x and y values 
  xs <- c()                                #and returns the borders of the 
  ys <- c()                                #Voronoi diagrams and where it should be centered
  for (i in seq_along(x_vals)){
    for (j in seq_along(unlist(x_vals[i]))){
      xs <- append(xs,unlist(x_vals[i])[j])
      ys <- append(ys,unlist(y_vals[i])[j])
    }
  }
  return(c((max(xs)+min(xs))/2,(max(ys)+min(ys))/2,max(xs),min(xs),max(ys),min(ys)))
}

gen_map <- function(coords){                              # Use those coordinates to generate the 
  return(get_map(location=coords[1:2],maptype='roadmap')) # NY map
}

The function generate_voronoi will plot a random sample of points into a Voronoi diagram on our map. We used ggmap to obtain the map previously

generate_voronoi <- function(df){                        # Plot Voronoi diagram onto
  df <- df %>%                                           # NY Map
    distinct()
  nymap <- gen_map(coords)
  ggmap(nymap)+
    xlim(coords[4],coords[3])+ylim(coords[6],coords[5])+
    geom_path(data=df, aes(x=Lon,y=Lat), stat="voronoi", alpha=.6) +
    geom_point(data=df, aes(x=Lon,y=Lat), color="red", alpha=.6)
}

The function find_pt takes in many points using this data to determine with polygon they are in.

find_pt <- function(xvals,yvals,df){                       # Determine how many given
  ex <- data.frame()                                       # points are within certain 
  for (i in seq_along(xvals)){                             # Voronoi polygons
    frame=data.frame(I(list(xvals[i])),I(list(yvals[i])))
    names(frame) = c("x_vals","y_vals")
    ex <- rbind(ex, frame)
  }
  ex <- mutate(ex,count=0)
  for (i in seq_along(df$Lat)){
    for (j in seq_along(ex$x_vals)){
      if (point.in.polygon(df$Lon[i],df$Lat[i],unlist(ex$x_vals[j]),unlist(ex$y_vals[j])) > 0){
        ex$count[j] <- ex$count[j] + 1
      }
    }
  }
  return(ex)
}

Finally, we use that data using create_choro to create a choropleth map of our data

create_choro <- function(df){                  # Put data from find_pt() to
  ids <- seq_along(df$count)                   # a choropleth map
  id=c()
  value=c()
  for (i in seq_along(df$count)){
    id <- append(id,rep(ids[i], each=length(unlist(df$x_vals[i]))))
    value <- append(value,rep(df$count[i], each=length(unlist(df$x_vals[i]))))
  }
  xs <- c()
  ys <- c()
  for (i in seq_along(df$count)){
    for (j in seq_along(unlist(df$x_vals[i]))){
      xs <- append(xs,unlist(df$x_vals[i])[j])
      ys <- append(ys,unlist(df$y_vals[i])[j])
    }
  }
  data <- data.frame(group=id,Count=value,x=xs,y=ys)
  nymap <- gen_map(coords)
  p <- ggmap(nymap) + 
    xlim(coords[4],coords[3])+ylim(coords[6],coords[5])+
    geom_polygon(data, mapping=aes(x=x,y=y,fill = Count, group = group), alpha=.70)
  return(p)
}

We took the past 6 months of Uber data from New York City, and we will create diagrams for each of them.

df <- read_csv(here::here('Uber-dataset','uber-raw-data-Sep21.csv'))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Date/Time` = col_character(),
##   Lat = col_double(),
##   Lon = col_double(),
##   Base = col_character()
## )
data <- slice_sample(df, n=10)
tiles <- deldir(data$Lon, data$Lat) %>% 
    tile.list
coords <- get_coords(get_poly_x(data), get_poly_y(data))
generate_voronoi(data)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.7777,-73.9397&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

get_area(data)
##  [1] 0.000360 0.000167 0.008988 0.002892 0.007003 0.001268 0.000150 0.000960
##  [9] 0.000638 0.008571
find_pt(get_poly_x(data), get_poly_y(data), slice_sample(df, n=25)) %>% 
  create_choro()
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.7777,-73.9397&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

df <- read_csv(here::here('Uber-dataset','uber-raw-data-Oct21.csv'))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Date/Time` = col_character(),
##   Lat = col_double(),
##   Lon = col_double(),
##   Base = col_character()
## )
data <- slice_sample(df, n=10)
tiles <- deldir(data$Lon, data$Lat) %>% 
    tile.list
coords <- get_coords(get_poly_x(data), get_poly_y(data))
generate_voronoi(data)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.74835,-73.9816&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

get_area(data)
##  [1] 0.000744 0.000717 0.000324 0.001057 0.000846 0.000224 0.000661 0.000685
##  [9] 0.000283 0.000258
find_pt(get_poly_x(data), get_poly_y(data), slice_sample(df, n=50)) %>% 
  create_choro()
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.74835,-73.9816&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

df <- read_csv(here::here('Uber-dataset','uber-raw-data-Nov21.csv'))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Date/Time` = col_character(),
##   Lat = col_double(),
##   Lon = col_double(),
##   Base = col_character()
## )
data <- slice_sample(df, n=10)
tiles <- deldir(data$Lon, data$Lat) %>% 
    tile.list
coords <- get_coords(get_poly_x(data), get_poly_y(data))
generate_voronoi(data)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.7467,-73.93685&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

get_area(data)
##  [1] 0.003767 0.000302 0.003062 0.000577 0.001906 0.001052 0.000314 0.000584
##  [9] 0.000415 0.000334
find_pt(get_poly_x(data), get_poly_y(data), slice_sample(df, n=5)) %>% 
  create_choro()
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.7467,-73.93685&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

df <- read_csv(here::here('Uber-dataset','uber-raw-data-Dec21.csv'))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Date/Time` = col_character(),
##   Lat = col_double(),
##   Lon = col_double(),
##   Base = col_character()
## )
data <- slice_sample(df, n=15)
tiles <- deldir(data$Lon, data$Lat) %>% 
    tile.list
coords <- get_coords(get_poly_x(data), get_poly_y(data))
generate_voronoi(data)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.734,-73.93585&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

get_area(data)
##  [1] 0.001306 0.000459 0.000637 0.001961 0.001895 0.001519 0.006726 0.000523
##  [9] 0.005876 0.002491 0.000385 0.001690 0.000229 0.003355 0.000926
find_pt(get_poly_x(data), get_poly_y(data), slice_sample(df, n=50)) %>% 
  create_choro()
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.734,-73.93585&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

df <- read_csv(here::here('Uber-dataset','uber-raw-data-Jan22.csv'))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Date/Time` = col_character(),
##   Lat = col_double(),
##   Lon = col_double(),
##   Base = col_character()
## )
data <- slice_sample(df, n=20)
tiles <- deldir(data$Lon, data$Lat) %>% 
    tile.list
coords <- get_coords(get_poly_x(data), get_poly_y(data))
generate_voronoi(data)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.7464,-73.8916&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

get_area(data)
##  [1] 0.001401 0.000362 0.000777 0.000598 0.000204 0.002454 0.000207 0.007277
##  [9] 0.000967 0.002783 0.000547 0.000606 0.005134 0.005484 0.010380 0.004358
## [17] 0.000029 0.006042 0.014375 0.000271
find_pt(get_poly_x(data), get_poly_y(data), slice_sample(df, n=75)) %>% 
  create_choro()
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.7464,-73.8916&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

df <- read_csv(here::here('Uber-dataset','uber-raw-data-Feb22.csv'))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `Date/Time` = col_character(),
##   Lat = col_double(),
##   Lon = col_double(),
##   Base = col_character()
## )
data <- slice_sample(df, n=75)
tiles <- deldir(data$Lon, data$Lat) %>% 
    tile.list
coords <- get_coords(get_poly_x(data), get_poly_y(data))
generate_voronoi(data)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.74685,-73.91035&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).

get_area(data)
##  [1] 0.000314 0.000335 0.001517 0.000017 0.000099 0.010611 0.015168 0.000026
##  [9] 0.000058 0.000363 0.000023 0.000294 0.002682 0.000024 0.000040 0.000654
## [17] 0.000038 0.000098 0.003541 0.003737 0.002586 0.000043 0.000078 0.000029
## [25] 0.000074 0.002390 0.000027 0.000580 0.000326 0.000307 0.000288 0.000063
## [33] 0.000023 0.004555 0.000038 0.000642 0.000020 0.000151 0.000112 0.001274
## [41] 0.000153 0.000088 0.000039 0.000533 0.000005 0.000122 0.000046 0.000049
## [49] 0.000033 0.000027 0.000148 0.000756 0.000181 0.000071 0.000056 0.000020
## [57] 0.000229 0.000023 0.000070 0.000011 0.000140 0.000164 0.000785 0.003246
## [65] 0.000033 0.000048 0.000030 0.000016 0.000037 0.002220 0.002391 0.000201
## [73] 0.000052 0.000052 0.008231
find_pt(get_poly_x(data), get_poly_y(data), slice_sample(df, n=1000)) %>% 
  create_choro()
## Source : https://maps.googleapis.com/maps/api/staticmap?center=40.74685,-73.91035&zoom=10&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx-nN5o2Fdo2txkNnbiUHbA
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 1 rows containing missing values (geom_rect).